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O Abstract 

> We numerically construct a family of five-dimensional black holes exhibiting a line of 

first-order phase transitions terminating at a critical point at finite chemical potential and 
temperature. These black holes are constructed so that the equation of state and baryon 
susceptibilities approximately match QCD lattice data at vanishing chemical potential. The 
critical endpoint in the particular model we consider has temperature 143 MeV and chemical 
potential 783 MeV. Critical exponents are calculated, with results that are consistent with 
mean-field scaling relations. 
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1 Introduction 



At zero chemical potential \x for baryon number, quantum chromodynamics (QCD) appears 
to have a smooth but rapid crossover at a temperature T c whose value is within about 
10% of 175 MeV. It is believed that this crossover sharpens into a line of first order phase 
transitions at finite /i. The position of the critical point that terminates this line is of 
considerable experimental interest, but it is hard to determine theoretically due to being 
in a region of strong coupling, and also because lattice techniques are not well adapted to 
finite real \x. A theory review can be found in [TJ. Recent lattice results can be found in for 
example [21 El H] . The aims of the recently initiated beam-energy scan at RHIC are laid out 
in [5] and the fixed-target CBM project at FAIR is discussed in [6]. 

It was shown in [TJ [8] that simple gravitational theories in five dimensions are capable 
of producing black holes which approximately reproduce the equation of state of QCD, 
including the crossover, at vanishing chemical potential; see also (9j EU [HJ H2J EE] and the 
review [TJ]. The gravitational theories include just two fields: the spacetime metric, and a 
real scalar field whose profile breaks conformal invariance and can be understood roughly 
as the running coupling of QCD. A natural generalization of such models is to include a 
chemical potential for baryon number as well. Adding a single additional field, a U(l) gauge 
field dual to the baryon number current, one may generate a chemical potential by turning 
on an appropriate electric field in the black hole geometry. 

In this paper, we study the coupled metric-scalar-gauge field system and numerically 
obtain solutions for charged black holes filling out the T-\x phase diagram of the field theory 
dual. The minimal Lagrangian for these theories contains some freedom, encoded in the 
choice of scalar potential and gauge kinetic function. We elect to fix both these functions 
by matching to lattice results for QCD at zero chemical potential. The scalar potential is 
fixed by demanding a QCD-like equation of state that captures the chiral symmetry breaking 
crossover, as in [7J|8]. We show how the gauge kinetic function can be determined by similarly 
matching quark susceptibilities, in principle removing all freedom from the construction. 

We then investigate how black holes behave at finite chemical potential. We find that 
just as is expected for QCD, at finite /i the crossover turns into a line of true first-order 
phase transitions ending in a critical point. We locate the first-order line by looking for 
characteristic thermodynamically unstable solutions, and identify the critical point as the 
end of this line; the location of the critical point is at physically reasonable values of T and 
/i. 

We then turn to a study of the critical exponents of this point. We find a set of exponents 
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that are nontrivially self-consistent due to satisfying two scaling relations. Thus our black 
holes built from just three fields reproduce a realistic phase diagram near the critical endpoint 
for a QCD-like theory. In QCD, the critical point is expected to lie in the same universality 
class as the 3D Ising model and the fluid liquid/gas transition. The critical exponents 
we obtain are consistent with mean-field scaling. This is reasonable since our black hole 
constructions are classical, corresponding to an implicit large N limit on the field theory 
side that suppresses quantum corrections. Further realism lies, presumably, in the inclusion 
of 1/N corrections. 

The organization of the rest of this paper is as follows. In section [2] we describe our 
gravity theory and summarize our results for the location of the critical point in the T- 
H plane and the values of its critical exponents. In section [3] we provide a self-contained 
summary of the aspects of thermodynamics which we will require in the rest of the paper, 
as well as brief remarks on the phase structure of QCD. Experts will have no reason to read 
this section, which contains no new results. In section [4] we analyze the equations of motion 
following from our gravity theory, explain how to extract thermodynamic quantities from the 
black hole solutions, and summarize our numerical strategy. In section [5] we will explain how 
the gauge kinetic function can be chosen to match lattice data for the baryon susceptibility 
at = 0. In section [6] we describe locating the critical point on the phase diagram, and 
in section [7] we analyze its properties, calculating the critical exponents and finding them 
consistent with mean-field scaling. In section [8] we compare our result for the location of the 
critical point to others in the literature, and conclude with some discussion. 



2 Gravity Theory and Summary of Results 



Our model falls in a class of five- dimensional gravitational theories including a real scalar 4 
and an abelian gauge field along with the spacetime metric, defined by the Lagrangian 
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(1) 



where we use mostly plus conventions. With energy dimensions assigned so that [k] = —3/2, 
[9 in] — = [0] =0, one finds that Q is almostQthe most general action using this field 
content one can have with at most two derivatives. Arbitrary functions of multiplying the 



1 The exception is that one could add a Chern-Simons term A A F A F, but this term has no effect on the 
classical equations of motion, which are what we will study; in addition it vanishes for the solutions we are 
going to consider, because they have only electric charge, not magnetic. Thus we neglect it. 
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Einstein-Hilbert and scalar kinetic terms can be removed by conformal transformations of 
the metric and reparametrizations of </>, respectively. 

The black hole geometries consist of metrics taking the form 

p 2B(r) 

ds 2 = e 2A(r) \_U r)dt 2 + dS 21 + f_ rfr 2 ( 2 ) 

L J h[r) 

along with an ansatz for the scalar field and electrostatic potential depending only on the 
radial coordinate r: 

= 0(r), A^dx" = <$>(r)dt . (3) 

The coordinates (t,x) cover Minkowski space, R 3 ' , while the radial coordinate r represents 
the holographic direction. 

As explained in [7] in the case of vanishing gauge field, a choice of the scalar potential 
V{4>) can be translated into a dependence of the entropy on temperature T. ([7J chooses 
to work equivalently with the speed of sound (? s = d log T/d log s.) In fact, if a desired 
dependence s{T) is specified, then — within certain limits — one can find the V((p) that leads 
to it. A reasonable fit, not too far from T c , to lattice results for s(T), is achieved with the 
simple choic^] [8J 

V(d>) = - 12CQSh ? + ^ 2 with 7 = 0.606 and b = 2.057 , (4) 

and L a constant related to the number of degrees of freedom. 

The black holes describing matter at finite chemical potential include a nonzero gauge 
field as well. This introduces the problem of specifying the gauge kinetic function f(4>). We 
can always use the freedom to rescale to set /(0) = 1. The matching of the speed of sound 
described in the last paragraph is completely insensitive to the choice of f(4>). However, /(</>) 
can be fixed if one knows the baryon number susceptibility at /x = 0. This susceptibility 
is in fact fairly well known from the lattice [2]. In this paper, we will not be systematic in 
finding V((f>) and /(0) through a fit. Instead, we will focus on the above choice of V((f>) and 
a similarly simple form for f(<f)), namely 



sech [f(0 - 2)] 
sech 



f(4>) = " , (5) 



5 

which as we will discuss in section [5] leads to susceptibilities in good agreement with lattice 
2 The constant 7 in J4| is unrelated to the critical exponent which will appear later in the paper. 



3 



results. 

It is probably impossible to find a string theory construction that leads precisely to the 
potential Q and gauge kinetic function ([5]). Thus we cannot claim that the theory (JIJ is 
dual to a specific known field theory. However, string theory constructions do typically lead 
to potentials which include sums of exponentials of canonically normalized scalars. Thus 
these functions are at least in the ballpark of expressions that can be derived from string 
theory, and it is reasonable to place the dual to our model in the broad class of strongly 
coupled, large- N gauge theories. 

Having fit V(4>) and f(4>) to lattice quantities at fi = 0, we are able to use black hole 
constructions to extrapolate outward into the T-fi plane, where we indeed find a critical 
endpoint. Through methods explained in sections [5] and |6j we estimate the location of this 
critical point to be 

T c = 143 MeV // c = 783 MeV . (6) 

We will compare this result to other estimates in the literature in the conclusions. 

Because we have not made a systematic study of the forms of V((p) and f(4>) that ap- 
proximately match lattice data at /i = 0, we are not in a position to provide theoretical error 
bars for the result ^ . It is best to view this result as a proof of principle that you can get 
a critical endpoint in the T-/i plane using AdS/CFT methods, and that the values ^ are 
within the theoretical error bars. It is also noteworthy that we ignore fluctuations in our 
analysis: the black holes we construct are fixed, classical geometries. This means that we 
are not capturing all the physics that is expected to go into the critical endpoint. 

Analyzing the thermodynamics near the critical point and performing linear regression 
fits to the data, we obtain results for four critical exponents, 

« = (), /9w 0.482, 7 ^ 0.942, 5^ 3.035. (7) 

These results are, as we shall discuss, non-trivially consistent with scaling relations, and 
consistent also with the mean field exponents a = 0, (3 = 1/2, 7 = 1, 5 = 3. 

3 Thermodynamics with a finite chemical potential 

Before discussing the solution of the equations of motion in the gravity system ([!]) and 
the exploration of the phase diagram, in this section we review a few essential aspects of 
thermodynamics and critical phenomena. 
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3.1 Thermodynamics of a fluid 

A fluid is characterized by the extensive quantities entropy S, volume V and particle number 
(or net charge) N, and their conjugate intensive variables temperature T, pressure p and 
chemical potential p. The internal energy U = U(S, V, N) depends on the extensive variables, 
which can be thought of as characterizing the system itself, with small changes described by 
the first law, 

dU = TdS - pdV + pdN . (8) 

The intensive variables, sometimes called "fields," can be thought of as properties imposed 
on the system by contact with a reservoir. 

For systems like the quark-gluon plasma, we are not interested in a fixed volume, but 
instead in volume densities for the extensive quantities. Define the energy density e, entropy 
density s and number density p, 

e = U/V, s = S/V, p = N/V. (9) 

Then, using the thermodynamic relation, 

U = TS-pV + pN, (10) 

one can show that the first law of thermodynamics ^ rewritten in terms of densities becomes 

de = Tds + \idp , (11) 

naturally reducing the system to a two-variable problem. It is useful to define the corre- 
sponding free energy density f(T, p) depending on the field variables T and p using the usual 
Legendre transformation]^] 

f(T,p)=e-sT-pp, (12) 

which obeys 

df = —sdT — pdp . (13) 



Furthermore, it is easy to see the thermodynamic relation (10) implies that the pressure 
reappears in the analysis as just minus the free energy, 

p = -f. (14) 
3 Thc free energy density should not be confused with the function f(cj>) in the gravity Lagrangian. 
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The phase diagram is the plot of "field" variables T and /i. At each point on the phase 
diagram, a physical phase corresponds to values of the extensive variables (or densities of 
extensive variables in our case, s and p) which extremize the free energy. In general, it is 
possible for more than one extremum to exist at a given point on the diagram, corresponding 
to the existence of multiple phases. The preferred phase is the one minimizing the free 
energy /; this is the condition of global stability and ultimately stems from the second law 
of thermodynamics. 

In addition to global stability, one must consider local stability, which is characterized by 
stability under small fluctuations. This is equivalent to the statement of positive-definiteness 
of the matrix of susceptibilities: 





/ ds 


ds\ 


)- 


dT 


dp 




\ dp 


dp 




\dT 


dp/ 



(15) 



where all derivatives of T or \x are taken with the other fixed. We note that the upper-left 
diagonal element is related to the specific heat at constant chemical potential C M : 

while the lower-right diagonal quantity is related to the isothermal compressibility, 

1 (dp\ 



^j\f,) T (17) 



In QCD, one usually considers in place of kt the quark susceptibility, 

dp\ fd 2 f 



^{r») T -{^)r p2KT - (18) 

The statement that the matrix of susceptibilities is positive definite is equivalent to requiring 

C P >0, X2>0, (19) 
where C p is the specific heat at constant volume, 



dT 2 (d 2 f/dfi 2 ) 



(20) 
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In the last expression all T and \i derivatives are taken with the other fixed. Note that the 
Jacobian of the susceptibility matrix is simply 

detS=^x 2 C p . (21) 

Generally speaking, local stability (i.e. positive definite S) is a requirement in order for a 
configuration to be considered a well-defined phase of the system. 

When two phases have equal free energies at a given (T, /a) , a first-order phase transition 
occurs at that point; typically the locus of first-order transitions has codimension one and 
so describes a line. On one side of the line one phase is favored, while on the other the other 
is favored; at the first-order line, a discontinuity exists in the densities, As and Ap, as the 
system jumps from one phase to another. The discontinuity in the entropy gives rise to the 
latent heat L = TAS. 

A first-order line may terminate on a second-order point, also called a critical point. 
Here the two distinct phases merge into one; consequently the discontinuities in the densities 
approach zero as one moves along the first-order line to the critical point. However, although 
these first derivatives of the free energy become well-behaved, the second derivatives, that 
is the susceptibilities, may diverge at the critical point. The nature of these divergences 
characterize some of the critical exponents associated to the critical point, as we will describe 
momentarily. Critical points in different systems with different variables may share the 
same critical exponents, a phenomenon called "universality;" the systems are said to lie in 
the same universality class. Beyond the critical endpoint on the phase diagram there is no 
discontinuous behavior, but the densities may change rapidly along the extension of what 
would have been the first-order line. This behavior is called a crossover. 

3.2 Phase diagram for QCD 

A great deal of work has gone into predicting the phase structure of QCD, and it is believed 
to be quite rich; for reviews, see [T5l [T6l [TT] . At small values of the baryon chemical potential, 
the QCD phase diagram is dominated by the chiral symmetry breaking transition, and this 
will be our focus. 

When all quarks are assumed to be massless, chiral symmetry is an exact symmetry of 
the QCD Lagrangian, and the broken symmetry phase at low T and p and the restored 
symmetry phase at high T and/or p are distinct and must be separated by a line of true 
phase transitions. Near the //-axis, the transition is expected to be first-order. Near the 
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Figure 1: The expected phase diagram of QCD. The line ending in a star is the first-order 
chiral transition and its critical endpoint, which we focus on. Below is the nuclear matter 
transition. At lower right are color superconducting phases, color-flavor locked and otherwise. 

T-axis, the order of the transition depends on the number of massless quarks. For two 
massless quarks, the transition on the T-axis is second order and in the universality class of 
the 0(4) model; this transition is expected to be the end of a line of second-order transitions 
extending into the T-fj, plane and meeting the first-order line rising from the /z-axis at a 
tricritical point. For three massless quarks, on the other hand, the transition on the T-axis 
is expected to be first-order. 

In the real world, quarks are massive and chiral symmetry is not an exact symmetry of 
QCD. On the T-axis, the transition is known from lattice studies not to be a sharp transition 
but instead a crossover. It is widely expected that at sufficiently large chemical potential /i 
the first-order line returns; it then terminates at a critical endpoint at some (T c ,/i c ). This is 
displayed in figure [TJ 

The critical endpoint is an object of substantial interest and speculation. It is difficult to 
explore it theoretically, as the theory is strongly coupled and lattice calculations are difficult 
at finite /i. A number of models have been constructed to analyze its properties. It is 
expected to lie in the universality class of the 3D Ising model, like the standard liquid/gas 
transition of fluids. It is anticipated that depending on its location on the phase diagram, 
future heavy ion experiments such as those at RHIC, LHC or FAIR may produce a quark- 
gluon plasma lying close to the critical point at freeze-out, which could lead to information 
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about its properties (see for example [51 El HE]-) 

Other phases of QCD are anticipated to exist, in particular regions at large p character- 
ized by color superconductivity. We will have little to say about these phases in this paper 



other than a brief speculation in section 6.2 



3.3 Critical behavior 

Near the critical point various first and second derivatives of the free energy go to zero or 
diverge as power laws, and it is the "critical exponents" associated with these power laws 
that are universal — meaning they may take the same values from one physical system to 
another, even among systems with quite different microscopic properties. Describing them 
is at the heart of the study of critical phenomena. 

We will calculate four standard thermodynamic critical exponents a, /3, 7 and In 
calculating the exponents, it is vital to specify whether one is approaching the critical point 
along the axis defined by the first order line, or by another direction. The exponent a is 
defined by the power law behavior of the specific heat at constant p as the critical point is 
approached along the axis defined by the first order line: 

C p ~ \T — T c \~ a , along first order axis . (22) 

The exponent /3 comes from the discontinuity of p across the first-order line. Ap is finite at 
a generic point on the first-order line, and goes to zero as one approaches the critical point 
along the line: 

Ap ~ (T c — Ty , along first order line . (23) 

The exponent 7 is analogous to a, but instead of C p , it is xi that is tracked along the 
first-order axis: 

X2 ~ \T — T c |~ 7 , along first order axis . (24) 

Finally, 5 is defined at the critical isotherm T = T c by the relation between p — p c and p — p c : 

p- p c ~\p- p c \ 1/s , for T = T C . (25) 

The same power law will manifest for any approach not parallel to the first-order line. The 
paths through the phase diagram associated to the four exponents are summarized in figure [2] 

Two other commonly-used critical exponents, v and 77, require knowledge of the spatial distribution of 
correlation functions and will not be calculated here. 
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Figure 2: A cartoon of the first-order line terminating at the critical point (star) with the 
directions of approach of the various critical exponents indicated. 

The four thermodynamic exponents are not all independent; in general they obey so- 
called scaling relations, which follow from the scaling behavior of the free energy at the 
critical point, determining two exponents in terms of the other two. One has 

a + 2/3 + 7 = 2 , 

(26) 

a + (3{l + 5) = 2. 

Different critical exponents are characteristic of distinct universality classes. Calculations 
in Landau-Ginzburg, or mean-field, theory capture the tree-level values of the critical ex- 
ponents; in general this neglects quantum corrections, which can be captured by the more 
sophisticated techniques of the renormalization group. The critical point of QCD is expected 
to lie in the universality class of the 3D Ising model, as does the standard liquid/gas transi- 
tion. The results from mean field (van der Waals) theory, the full quantum 3D Ising model, 
and experiments in non-QCD fluids are summarized in the table [PD] : 





Mean field 


3D Ising 


Experiment 


a 





0.110(5) 


0.110 - 0.116 


& 


1/2 


0.325 ± 0.0015 


0.316 - 0.327 


7 


1 


1.2405 ± 0.0015 


1.23 - 1.25 


5 


3 


4.82(4) 


4.6 - 4.9 



These are the results we will compare our holographic system to. 
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4 Black hole solutions 



We now turn to an analysis of the equations of motion for the gravity system. From the 
action one can derive four second order equations of motion and a zero-energy constraint. 
The second order equations (simplified slightly using the zero-energy constraint) are 

A" - A'B' + V 2 = 
6 

h" + (AA' - B')ti - e- 2A f((j))& 2 = 

$" + (2A' - B')& + = (27) 

where 



Vafa r) = V{<P) - y^m®' 2 . (28) 
The zero-energy constraint is 

h(2AA' 2 - (j)' 2 ) + QA'ti + 2e 2B V(<P) + e- 2A /(0)$' 2 = . (29) 

The equation of motion for $ can be integrated once to show the conservation of the Gauss 
charge Qg for the U(l) gauge field: 

^=0 where Qg = f((j))e 2A ~ B & . (30) 
dr 

One other conserved quantity can be guessed from scaling symmetries, as in [7]: 

^ = where Q N = e 2A ' B [e 2A h' - f (<{>)$&'] . (31) 
dr 

4.1 Near- horizon asymptotics 

Let's assume that h has a simple zero at r# and that it has no additional zeroes between r# 
and the boundary. Then r# is the location of a regular black hole horizon. A series solution 



to the equations (27) and (29) can be developed simply by expanding 



X(r) =X + X l (r~ r H ) + X 2 (r - r H ) 2 + . . . , (32) 
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where X is any of A, B, h, $, and <fi. B(r) may be fixed to be anything by a choice of 
the coordinate r, so all the B n are arbitrary. All but finitely many of the other coefficients, 
however, are determined in terms of the first few. To be more precise: ho = by assumption; 
A = can be arranged by rescaling t and x by a common factor; hi = 1/L can be arranged 
by rescaling only t; $ = is a choice one must make in order for Qdt to be well-defined at 
the horizon; and all other coefficients are determined once one chooses 0o and In other 



words, the solutions to (27) and (29) may be parametrized by (0o, 3>i)- It helps our intuition 
to recall that 0o is the value of the scalar field at the horizon, while $1 is essentially the 
electric field in the radial direction, also evaluated at the horizon. Given the assumptions 
just stated, it is easy to show that 

Qg = e- B V(0o)$i 

1 _ (33) 
Qn = jf B * ■ 

It does not seem to be practical to find solutions of the equations of motion through high- 
order series expansions, because the expressions for high-order coefficients quickly become 
quite complicated. In practice we stopped at fourth order. The resulting expansions are 
suitable for providing initial values for numerical integration of the differential equations 



(27) at a radius slightly outside the horizon. 



4.2 Far region asymptotics 

To discuss asymptotic behavior far from the horizon it helps to pick a gauge, so we fix B = 0. 
We can write the potential V((p) as 

V(<f>) = -^ + lmy + 0(<f> 3 ), (34) 

demonstrating L is the radius of curvature of the asymptotic AdS^ geometry, and we can 
define A^, the ultraviolet dimension of the operator dual to 0, according to 

m\l? = A (A^ - 4) . (35) 

Following [S] , we have assumed that this operator is essentially tr F 2 and that the ultraviolet 
limit is to be matched to QCD at a scale significantly above T c but not parametrically large. 
Thus A^ should be only slightly less than 4, and in our potential A^ w 3.93. 



12 



One can then straightforwardly show that 



where 



and 



A(r) = a(r) + A^e~ 2va ^ + . . . 

h(r) = h l Q ai + ^ ar e" 4a(r) + h% 2v e-^ 2v)a(r) + ... 

$( r ) = $f ar + $f?r e -2 a (r) + $ far^ e -(2+^)a(r) + _ _ _ (35) 

0(r) = <p A e~ ua{r) {l + a v e- va W + a 2v e~ 2ua ^ + . . .) 
+ <j) B e- A * a{r) + ... 



a (r) = A^y + (37) 



i/ = 4-A*. (38) 



The zero-energy constraint (29) implies 

1 



A! ' 1 = 7W 



Given (36), it is straightforward to show that the conserved charges in terms of the far-region 



quantities evaluate to 

Qg = -jA^f 

o L (40) 

Q N = A^\{-2hf r + ■ 
Lj 

Thus A f ^[, h f £ r , and $2 ar can be determined in terms of h^ r and $Q ar once Qg and Qn are 
known. It is notable that in the absence of a scalar (or if for some reason — > at the 
boundary faster than e~ 2a ^) then the next correction to h(r) after ft,| are_ ^ is /ig ar e -6Q ( r \ 
and the equation of motion for h can be used to show that 

hi ,r = 1(^)2 (41) 

However, for the solutions we will study numerically, the scalar doesn't vanish fast enough 
for the h}§ r e~ &a ^ term to be interesting. 



Evidently, the expansions (36) are qualitatively more intricate than the Taylor expansions 
(32) around the horizon because the powers of e~ a ^ are not (for practical purposes) com- 
mensurate, owing to v m 0.07 not being the ratio of small integers. The fact that v <C 1 also 
leads to some difficulties in finding robust numerical solutions to the equations of motion. 
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We will discuss these issues at greater length in subsection |4.4 



The expansion of <p{r) in (36) is split into the part dual to a deformation (proportional 
to 4>a) and the part dual to an expectation value (proportional to 0b). Each solution carries 
a series of corrections, which we have shown only for the solution proportional to (f>A- The 
correction term a„e- va ^ is present only when V'"(0) 7^ 0, so for even potentials like Q the 
leading correction is a2 U e~ 2ua ( r \ For small enough v, this correction, and even higher cor- 
rections proportional to (J)a, dominate over the (j)Be~ A < t ' a ( r ' ) term. Thus the terms not shown 
explicitly in the expansion for 0(r) are subleading either to 



to both. In all the other expansions in (36), the omitted terms are all subleading to the 
terms shown explicitly. 

4.3 Thermodynamic quantities 

The solutions we are interested in have 4>a 7^ 0, because they are to be understood as renor- 
malization group flows triggered by deformation of a very slightly relevant operator. This is 
not exactly how QCD works, but it is sufficiently close to be an interesting approximation. 
In order to compare solutions meaningfully, one should ideally arrange for <pA always to be 
the same. This can be accomplished through a coordinate transformation, provided 0(r) 
always has the same sign at the boundary]^] More specifically: suppose we obtain a solution 
numerically which has some positive value of <pA- Then we wish to perform a coordinate 
transformation on this solution to bring it into the form 

d ~2 = e *Mf)(-h(f) d p + d &) + ifL 

h(r) (42) 
A^dx^ = $(r)dt 4> = 0(f) 

where 

A(r) = y + 0(e- 2vf/L ) 

h(f) = 1 + hfe- Af ' L + (e-^ +2 ^ L ) (43) 
$(f) = $f ar + ^ T e - 2f/L + 0(e- {2+u)f/L ) 
0(f) = e ~ vf/L + 0(e~ 2vf/L ) 



5 If <fi(r) becomes negative at the boundary, the following expressions for far-zone coefficients and thermo- 
dynamic quantities can still be used if \4>a\ is substituted for <pA- 
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Setting ds 2 = ds 2 , A^dx^ = A^dx^, and <f>(r) = 0(f), one finds immediately that 



x 



x 



(44) 



and 



a(r) - log^Y") 

i(f) 
h{r) 



A^\- + A Q 



far 



log( 



A(r 
1 

/,far 
n 



h(r) 



(45) 



$(r) 



which implies 



far 



$ far 



l. 1 / 1 ' /Tfar 
A V n 



far 



/if 



$ 2 ar 



(46) 



/if 



/4/^far ' 



'A n 



Intensive thermodynamic quantities can now be readily extracted using these relations along 
with standard expressions in the (t, x, f) coordinate system and the assumptions stated 



following (32): 



T 



Air 




r=r H 



4?r l$H v 



a 



$ far 



(47) 



L I4]l v y/f¥' 

Densities of extensive thermodynamic quantities can likewise be computed: 



2tt 



2tt 1 



K 



P 



(48) 



K 



far 




Thus knowledge of the four asymptotic scaling parameters 4>a, h^, $o ar and $2 determines 
the standard thermodynamic variables T, /1, s and p. Subleading parameters in the field 
expansions will depend on these in general. For example, using the constancy of the Noether 
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charge (|3 1 h and its asymptotic expressions (]33|) and (]40|), for ft,| ar one can show 



hf 1 



h 



far 



K 2 L 



{sT + up) 



k 2 L 



e + p) 



(49) 



where in the last step we used the thermodynamic relation (10) to express the result in terms 



of the sum of the pressure and energy density. Note that this is the only combination of e 
and p we have access to from these calculations. The pressure by itself is equivalent to the 



free energy density (14), which to calculate we would need to evaluate the full renormalized 



action including counterterms to cancel divergences. It is possible to get the results we're 
interested in — in particular the position of the critical point and the values of its critical 
exponents — with just T, p, s, and p. 

We can also use the Gauss charge to relate a certain combination of the near-horizon 
parameters (<fto, $i) to the asymptotic parameters, and thus the thermodynamics. One finds 
that the Gauss charge is proportional to the inverse of the entropy per baryon: 



Q G = /(0 o )$i 



47T p 

Is 



(50) 



This is the only analytic relation between the initial conditions at the horizon and the 
thermodynamic parameters. 



4.4 Numerical strategy 

It is straightforward in principle to obtain a numerical black hole solution by integrating the 



second-order equations of motion (27) starting at some point slightly outside the horizon with 
the functions and their derivatives initialized from the horizon series expansions described in 
section 4.1, with initial conditions (0o> $i)- Then a fit can be performed of the numerically 



known functions A(r), h(r), ^(r), and 0(r) to the asymptotic forms (36) to extract the 
quantities h^, $Q ar , $2 ar , and <$a in terms of which T, p, s, and p can be determined. Thus 
for each input value of (0o? ^i) ? we obtain a black hole characterized by thermodynamic 
quantities (T, p, s, p). Certain values of (</>o, $i) may lead to a solution that does not converge 
to an asymptotically- AdS^ solution. Typically, these spacetimes are singular. They are not 
of the class we are interested in, so they are discarded. 

Numerical integrations can be made vastly more efficient by noting that h(r) and $(r) 
converge much faster to their asymptotic values than 0(r) and A'(r). A good strategy, then, 
is to figure out the value r = r* beyond which the non-constant corrections to h(r) and $(r) 
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have no more influence on the equations of motion for A and than round-off errors do; 
then join a solution of the full equations of motion from a point just outside the horizon to 
r = r t to a solution to simplified equations of motion, obtained by replacing h by and 
$ by $o ar , from r = r t to a value of r large enough to reliably compute <pA- As discussed 
following (40), $2 ar can be determined once h l Q V and $Q ar are known. Thus in order to extract 



T, p, s, and p using (47)-(48), the only quantities one needs from numerics are h^, $q, and 



^a- We implemented the strategy described here in Mathematica, where the basic ODE's 



(27) are solved using NDSolve. 



5 Quark susceptibility at zero chemical potential 

Because lattice calculations at finite chemical potential are problematic, it has been difficult 
to make precise predictions for the behavior of QCD off the T-axis. However, at p = 0, lattice 
studies have been carried out extensively. The potential V(4>) from [TJ |8] was engineered to 
reproduce the equation of state s(T) known from lattice simulations. 

We would like to also constrain the gauge kinetic function /(</>) using known lattice results 
at = 0. The extrapolation to finite p is then completely determined by known physics 
at p — 0, and represents the unique prediction for the phase diagram of the large- N gauge 
theory defined to emulate the thermodynamics of QCD on the T-axis. 

The gravity calculation of s(T) at p — is completely insensitive to /(</>), since the 



gauge field is zero in these solutions. Instead we may examine the quark susceptibility (18) 



at vanishing p as a function of temperature, as this has also been calculated extensively on 



the lattice and as we will see, depends on the choice of f(4>). In section 5.1 we find a gravity 
formula for the quark susceptibility at zero chemical potential, and in section 5.2 we use this 
to justify our choice of f(<fi). 

5.1 A formula for quark susceptibility 

The black holes with p = have vanishing gauge field A^, and p = as well. To calculate the 



quark susceptibility ( 18 ), we make use of the key observation that the gauge field equation of 



motion is linear and homogeneous in $, while $ appears only quadratically in the remaining 



equations (27). We thus proceed by treating $ as a linear perturbation, solving the gauge 



field equation in the fixed background of the p = black hole, and then determine \2 by 
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noting that on the T-axis, its definition ( 18 



) becomes 



X2 (^ = 0) = lim 



(51) 



Moreover, in the linearized approximation, the overall normalization of $ is arbitrary as far 



as the equations of motion are concerned and will cancel out of (51), so we can set it to 
$i = 1/L. 



We can in fact obtain a formula for (51) that reduces to quantities only involving the 
metric and scalar, which are unchanged in the linearized approximation and thus can be 
taken from the solution for the background fi = black hole. Since it is common in the 
literature to plot \2 normalized by T 2 , which approaches a constant at large T, we will find 



a formula for % 2 = X2/T 2 . Using equations (47)- (48), we have 

p (4vr) 2 L 3 $!f r ^ ar 



X2(/4 = 0) 



/jT 2 



(52) 



and making use of the expression (40) for the Gauss charge Qg, we may simplify (52) to 

^ 2 -k 4 ^fan3/2 Qg 



(53) 



Now recall that $ — > at the horizon r = r#. As a result, 



POO /'DC 

$f-= / dr& = Q G dre- 2A f((f>)- 1 

Jr H Jr H 



(54) 



where in the second step we have employed the definition (30) of the Gauss charge. Plugging 



(54) into (53) results in 



8tt 2 L 4 



far \ 3/2 



k 2 f™dre- 2A f(<P)- 1 



(55) 



Note at this point that all explicit dependence on the gauge field $ has dropped out. One 
can illuminate this further by noting that 



128tt 4 L 3 



T 3 



far\3/2 



so that finally 



Xa(^ = 0) 



16tt 2 T 3 J^dre-^fi^)- 1 



(56) 



(57) 
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This expression may now be evaluated on /x = black holes directly, without having to solve 
the linearized $ equation explicitly. 



Our final expression (57) is suggestive because in lattice simulations, s/T 3 and %2 = 
X2/T 2 have qualitatively similar behavior as functions of temperature at fi — 0: both start 
near zero for low temperatures, then rapidly cross over to a large value in the region of T c , 



and asymptote to a finite value at large T. Hence from (57) we come to expect the realistic 
behavior of X2 will to some extent be inherited from the analogous behavior of s/T 3 . 

The effects of the integral in the denominator of (57) do play an important role, however. 
This integral introduces a dependence of the quark susceptibility on the function f(4>), which 
s/T 3 alone was insensitive to. Thus differences between the functional forms of X2 and s/T 3 
are due entirely to the effects of /(</>). 



5.2 Matching to lattice data at zero chemical potential 

We now discuss the matching of lattice data to black hole results at ji = and justify our 
choice ([5]) for /(</>). Since the precise field theory dual of our model is unknown, we will 
not try to translate the quantities k, L into field theory language. Instead, we will make 
the arbitrary choice k — L — 1 and parametrize our ignorance by allowing separate overall 
constant rescalings between the lattice quantities and the black hole quantities: 

I ■'J lattice A s [s]bh, [T] lattice At[T]bh, [p] lattice lattice A M [/i] B H- ( 58 ) 

As described previously, to compute the entropy density at \i = one doesn't need any 
information about /(</>) at all. In figure [3]A. we show how the entropy density compares 
between lattice and black holes based on the potential Q. For lattice data we used the right 
hand plot in Figure 3 of [2j^J with points from asqtad N T = 6 simulations from T = 150 MeV 
out to T = 382 MeV, and then points from p4 N T = 6 simulations out to T = 720 MeV. We 
determined T c 191 MeV as the temperature at which s/T 3 reaches 1/e of its largest value 
as obtained from the highest temperature data point Q For black hole data, we constructed 
black holes starting with our standard choice Q of scalar potential, and for <pQ ranging from 

6 We chose [2J to have a definite lattice result to compare to, but there is still disagreement in the literature; 
for another determination of the equation of state, see @J. We expect small changes to our model could 
accomodate variations in the lattice results. 

7 T C 191 MeV is somewhat larger than the value T c ss 175 MeV mentioned in section [TJ however it is in 
line with estimates of [2D]. Lower values for T c are favored, for example, in J2TJ [22] [23] . We do not aim here 
to probe the apparent discrepancy; instead we are largely opting for the higher values because all the lattice 
data we use directly is from [JJ. 
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Figure 3: The normalized entropy s/T 3 and quark susceptibility x 2 = X2/T 2 at p = 0, 
computed on the lattice and fit by black holes in the gravity theory defined by our choices 
of V((f)) and /(</>) (equations Q and (|5|). Lattice data is taken from [2]. 



1.5 to 7.5 in 20 steps, uniform on a log scale. 

Our conventions are for all lattice quantities to have units which are powers of MeV, 
while with k = L = 1, all black hole quantities are dimensionless. Thus A s and Xt have 
units which are also powers of MeV. We found a good fit between lattice data and black 
holes with 



X s = (121 MeV) ; 



A 



T 



252 MeV 



(59) 



Turning to the quark susceptibility, in figure |3p we show how susceptibilities computed 
starting from the choice ^ for f((j>) compare with lattice data. For lattice data we used 
the same reference as for the entropy [2], with light quark results from the left hand plot in 
Figure 5, and strange quark and total baryon number results from the two sides of Figure 6; 
we scaled the light quark and strange quark curves appropriately to asymptote to the same 
value as that of the baryon number, X2 = 1/3, at high temperatures. For black hole data 
we used the same black holes as in the entropy plot, and we employed (57) with L = 1 (and 



implicitly also k = 1 as before). We used the value of At in (59) to rescale the temperature 
axis, and we adjusted the overall scale of \i arbitrarily to optimize the fit to lattice over the 
range shown in the figure [3)3. The rescaling of the susceptibility is thus 



be 



2 BH 



dp 
dp 



A 



2 I lattice 



(60) 



BH 



Thus, knowing [x2]bh and [X2] lattice & t the same temperature tells us A M /A p . In order to find 
A p and X p separately, we must recall that the relation for the free energy (12) holds equally 
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in lattice units and in the black hole setup. Thus 



\ T \ S = X X = \ e = \ f . (61) 



Putting (60) and (61) together, we find 

[X2]BH 

which can be recast as 



A?, , , Xj'X^ 



7< 



[X2]lattice — — r~ ~[X2] lattice > (62) 



\2 \ \ LA .4 J lattice , 
A T ApA^ A 



A s [xs 



A^Wfr^- (63) 

AT [X2\ lattice 



Let us now describe how we arrived at the choice ^ for functional form for the gauge kinetic 
function. The value of /(</>) near the horizon is particularly important, because the factor 
e -2A j n integral (57) puts significant weight on the near- horizon region. In particular, 



if f(4>) is large at the horizon, the integral will be relatively small compared to when /(</>) 
is small at the horizon. Since \2 stays close to its high-temperature value down to a lower 
temperature than s/T 3 before plunging rapidly to small values, a reasonable conjecture is 
that as <p goes from to positive values, one needs f(<f>) first to increase as a function of 
(p, then to decrease rapidly. The functional form ^ was chosen with these desired features 
in mind, and also with the thought that asymptotically exponential behavior at large is 
typical of supergravity theories. 

Operationally, the way we determined A M was to use the correct Stefan-Boltzmann value 
X2 = 1/3 for baryon number as the value for [£2] lattice; an d to evaluate [x2]bh at T = 460 MeV. 
This is a reasonable approach because the lattice data converges quickly to the Stefan- 
Boltzmann value at high temperature. The result is 

Xp = 972 MeV , A p = (77 MeV) 3 . (64) 

Again it should be emphasized that our choice ^ of /(</>) is to a degree ad hoc, and it should 
be understood as providing a proof of principle that an approximate fit to \2 can go with a 
critical endpoint in the T-jj, plane based on AdS/CFT techniques. 
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6 Searching for the critical point 



Having settled on a functional form for V(<f) and /(0) by matching to lattice thermodynamics 
at /i = 0, our Lagrangian is now completely determined. We can next turn to numerically 
solving for a set of black holes to fill in the phase diagram. An expectation is that the 
crossover that takes place on the T-axis is sharpened into a first-order line lying out in the 
T-fi plane, and that this first-order line ends at a critical point somewhere in the vicinity of 
the crossover. Our first task therefore is to search for this critical point. 



6.1 Scanning the thermodynamics of black holes 

For a first pass at mapping out the thermodynamic behavior of black holes across the T-fi 



plane, we generated approximately 2500 numeric solutions to the equations (27) and (29) 



seeded by initial conditions near the horizon as described in section |4.1| Each solution is 
specified by the value of (cj) ,^i) that was used to generate the near-horizon asymptotics. 
We remind the reader that <f) is the value of the scalar field at the horizon, which we took 
always to be positive, and $i is essentially the electric field at the horizon pointing upward 
in the fifth dimension. We worked exclusively in the gauge B = 0. 

To choose a suitable range for O , we note first of all that the fits discussed in the previous 
section involved values of 0o no larger than 7.5. We went from O = 1 to 0o — 15 in order to 
obtain the best global picture of the thermodynamics that we could, but any features seen 
at <po significantly larger than 7.5 should be regarded with some degree of skepticism, since 
in principle one could adjust V(4>) and/or /(</>) for <ft > 7.5 to make any desired phenomenon 
occur in that regionj^] It is notable, however, that both V(<j)) and /(</>) are fairly featureless 
for 0^4, both being close to a simple exponential function over that domain. 

To choose a suitable range for $ 1; we demonstrate that there is an upper bound on 
possible $! values leading to an asymptotically-AdS black hole. To see this, first note that 



the first equation of (27) with B = shows that A is concave down as a function of r. 
But it must be increasing at large r in order for the spacetime to be asymptotically AdS*>. 
Therefore A must be increasing at the horizon, which is to say A\ > 0. Using the zero-energy 



constraint (29), A\ can be re-expressed as 



A 1 = ~ [2Vfo) + /(to)*?] • (65) 



8 In fact, constraints on f(<p) come from a narrower range of <f>, extending only up to (f) = 5. Thus, 
baryon-specific physics is most reliably studied in our model at values of 0o no greater than 5. 
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Figure 4: Numerically generated black holes. Each dot represents a numerically generated 



solution. If the Jacobian J defined in (67) is positive for this solution, then the dot is red. 



If J < 0, it is green. The bold black circle is the critical endpoint. 



Because V(<f>o) < and /(0o) > 0, this puts an upper bound on $i: 



|$i| < $i 



' /(0o) 



(66) 



In practice we scanned black hole solutions from $x just slightly greater than up to 

Figure [4] shows the results of our numerical scan of the T-\i plane. We examined 61 values 
of 0o between 1 and 15, uniformly spaced on a log scale. For each value of 0o so obtained, 
we examined 41 evenly spaced values of $i/$i, m ax- A small fraction of the values so chosen 
failed to produce good black hole solutions, generally because A failed to be monotonically 
increasing, and are simply omitted from the plots. 



6.2 Locating the critical point 

To locate the critical point, we must think a little about what we expect to find in the vicinity 
of the first-order line. When there are competing phases in a thermodynamic system, only 
the one minimizing the free energy is the true ground state. However, there is no reason 
to think our black hole solution-generating method will discover only the true ground state 
solutions. Due to chiral symmetry not being exact and the presence of the crossover, there 
is no invariant distinction between the two sides of the first-order line that could correspond 
to a difference in topology or other invariant distinction between the phases on the gravity 
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Above T c At T c Below T c 




Figure 5: The baryon density p as a function of chemical potential p for several values of 
T near the critical point. For T > T c , p(p) is single- valued (left), while for T < T c it is 
multi- valued (right). At T = T c the slope is infinite (middle). 



side; the distinct phases will be continuously connected in the space of solutions. Since we 
are just solving the equations of motion, we expect to find all extrema of the free energy. 

In general, extrema of the free energy include not just locally stable minima, but also 
any thermodynamically unstable saddle points or maxima. Far from the first-order line on 
the T-p plane we expect only one solution to the equations of motion; in the vicinity of the 
first-order line, however, we expect to find p and s to be multivalued. Note that this will be 
true not only on top of the first-order line, but also merely near it, as the free-energetically- 
unfavored phase will persist for some distance on the phase diagram before ceasing to exist as 
a solution. The first-order line ends precisely at the (T c , p c ) where this multivalued behavior 
ceases; this is the critical point. 

Thus if we consider a constant-T slice of the phase diagram with T > T c and vary p, this 
isotherm will miss the first-order line and the functions p(p) and s(p) will be single- valued 
(although for T close to T c they will display crossover-type behavior). But for T < T c , the 
isotherm will intersect the first order line and we expect p(p) and s(p) to be multivalued 
near p c . The simplest behavior that still increases at both large and small p is an "S"- 
shape, and this is what we observe; see figure [5} For such behavior there are three solutions 
at a given p. Since the slope of the curve is just the quark susceptibility (18), we see 



that two of the solutions have X2 > and thus may be thermodynamically stable (19); 
these are the candidate phases. The middle solutions, however, have x% < and must be 
thermodynamically unstable]^] Precisely for T = T c the curve will cease to be multivalued, 
as the three solutions coalesce into one; the curve p(p) will have an infinite slope, indicating 



9 According to the correlated stability conjecture (CSC) [24] 25] . such black hole solutions will also have 
dynamical instabilities, corresponding to the black hole gaining total entropy by locally redistributing charge 
and energy subject to global conservation of these quantities. In the black hole literature this is known as 
the Gregory-Laflamme instability [551 HZ] ■ 
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a divergence in the quark susceptibility at the critical point. 

We will locate the critical point by looking for the thermodynamically unstable solutions 
that characterize the vicinity of the first-order line. We can identify the unstable solutions 



by calculating the Jacobian of the susceptibility matrix (15): 



J = detS = d{s,p)/d{T,fi). (67) 



For a thermodynamically stable black hole, equation (19) is satisfied and the Jacobian (21) is 
manifestly positive. If it flips sign to J < 0, we have necessarily found a thermodynamically 
unstable branch. Once we find the thermodynamically unstable black holes, we look to see 
whether they map to a narrow line-like region on the T-/i plane; the critical point is then 
the values (T c , /i c ) where this line ends. We should also be able to see the two stable phases 
mapping to the same locus on the phase diagram from elsewhere in (0 O > 

We can calculate the Jacobian J by finite differences. Since the black holes were scanned 
on a rectangular grid, we can label them with indices ij, where i determines the value of 0o 
and j determines the value of $i/$i, m ax- In order to compute J for the black hole labeled 
ij, we first computed 

jTfi _ ^ ^ I — Tij Tij + i — Tij 

^Pi+l,j Pi,j Pi,j+1 Pi,j y 

(68) 

,Pi+l,j — Pij Pi,j+1 ~ Pi,j / 



Then Jf? '/ ' jT^ is the finite difference approximation to the Jacobian J in (67). The results 
are shown in figure |4} 

It is clear from figure |4|\ that there is a region of unstable black holes stretching down to 
$i ~ 0.4$i >max . It is this region which we are most interested in, because when mapped to 
the T-fi plane it becomes a narrow region that ends in a cusp. Moreover, we do indeed find 
two other sets of black holes with J > mapped to the same locus, and hence we identify 
it as the first-order line. The point of this cusp is then the critical endpoint, which we show 
as a bold black circle. It occurs at the values 

(T c , n c ) w (143 MeV, 783 MeV) . (69) 



We have used the multipliers At and from (59) and (64) to express T and /x in units of 
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MeV. Points at the critical point come from initial conditions in the vicinity of 

(0o, $i/$i, max )~ (4.84, 0.40). (70) 

Note that the value for O is within the range probed by the /i = solutions described 
in section |5.2] — though not by much, if one goes by the values of <fi over which /(</>) is 
meaningfully constrained by lattice data. 

In summary, we have identified a candidate critical point and first-order line. As we study 
it in more detail in the next section, examining the behavior of densities and susceptibilities, 
this identification will be amply confirmed. 

Before moving on, let us consider the other thermodynamically unstable black holes found 
in our scan. The unstable black holes described so far are associated with a failure of the 
map (0o? ^i) ~~ >* (T,fi) to be invertible, which is to say a sign change in jj^. Only with 
such multiple covering can you jump abruptly from one solution to another at the same 
(T, n) but different (s,p): the sine qua non of first-order phase transitions. As one proceeds 
further to the right in the T-fi plane, one encounters a broader region of unstable black holes 
immediately above the multiply-covered region. This region is unstable due to a change of 
sign in J*?, meaning it is the map of (0o, ^i) — > (s, p) that is not invertible. This sign change 
causes black hole instabilities, presumably of the Gregory-Laflamme type (261 El], but no 
first-order line. Correspondingly, there are no stable black holes in this region of the phase 
diagram, 10 

The absence of stable black holes in our model at large fi (roughly, larger than fi = 
1100 MeV) and T not too big is actually a good thing. It is in approximately this region 
that one might reasonably expect color superconductivity and/or related phenomena to set 
in: see for example Figure 7 of [2H] and Figure 1 of [3] . Black holes based on the lagrangian 
([!]) are not likely to capture such phenomena. However, it is comforting to note that in cases 
where black hole superconductivity is understood (where the condensate breaks a U(l) gauge 
symmetry in the bulk), the superconducting instability competes against Gregory-Laflamme 
instabilities, and one generally must pass beyond minimal supergravity lagrangians to see 
the superconducting instabilities: see for example [29l [30] . Thus all our findings are at least 
qualitatively consistent with consensus expectations for the QCD phase diagram. 

10 It is interesting to note that J"? tends to change sign close to $i/$i, mM ~ 0.6 for a fairly wide range of 
other choices for f((f>) that we examined numerically. 
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7 Analysis of the Critical Point 



Having found the critical point, our final task is to determine its critical exponents. To 
achieve this, we first construct a large data set which densely populates the critical region. 
This collection of about 120,000 black holes is generically described by solutions with O £ 
[4.25,5.5] and $i/$i, max G [0.35,0.43]. 

Using these near-critical black holes, we can systematically study the approach of various 
thermodynamic quantities to criticality. Since the behavior of these quantities is typically 
expected to be power law in the vicinity of the second order point, it is natural to study 
them on log-log plots on which critical exponents are trivially related to the slope of the 
best fit to the data. In practice, we extract this slope by performing a linear regression via 
a least squares fit. All reported critical exponents in this section have been obtained in this 
way. 

7.1 First-order line and critical density 

Despite knowing the location (T c ,/i c ) of the critical point, it remains a challenge to identify 
the critical density p c \^\ Due to the infinite slope of the curve p(p) on the critical isotherm, a 
large number of black holes with very different values of p sit very close to the critical point 
(see the middle plot in figure [5]) . 

We can calculate the critical density in the process of determining the exponent /3 which 
measures the rate that the discontinuity Ap across the first-order line goes to zero as the 
critical point is approached: 

Ap ~ (T c — Ty, along first order line . (71) 

In order to determine the discontinuity in p for a given Tf < T c , we must identify the pj at 
which the true ground state of the system jumps from the lower branch to the upper branch; 
this is the pf for which the free energy is the same for the two branches. However we have 
not calculated the free energy, so we cannot identify pf in this way. Equivalently, one can 
use Maxwell's equal-area construction, which states that pf should be placed such that the 
closed regions bounded by the isotherm on either side of the p = pf line are of equal area. 

We chose a computationally easier procedure which is asymptotically equivalent to the 
equal-area law as one approaches the critical point. Namely, at a fixed temperature T = Tf, 

11 Remarks about p in this section apply equally well to s. 
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Figure 6: Cartoon of p(p) for an isotherm with Tf < T c , showing multivaluedness near the 
first-order line. At the location of the line p = pf the true minimum of the free energy jumps 
from the lower to the upper branch and p is discontinuous. 



we define /i< and p > to be the locations of the local minimum and maximum of the isotherm 
p(p), and we define pf to be the midpoint between them. This in turn determines the mixed- 
phase densities p< and p> for the point (T/,pf) along the first-order line. This procedure is 
illustrated in figure [6} For points we checked near the critical point, this procedure agrees 
with the equal-area rule to within a fraction of a percent. 

The critical density p c is then most easily obtained as the limit that both p < and p> 
approach as we near the critical point. The result is 

p c = 9.9022 . (72) 

Plotting Ap = p > — p< in a log-log plot with t=(T — T c )/T c , we obtain 

P w 0.482 . (73) 

In comparison, the exponent in the mean-field case is Pmf = 1/2, and so we have found a 
result very close to the mean-field value. 
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Figure 7: The discontinuity in the baryon density as the critical point is approached on a 
log-log plot. The slope of a best fit line through the data gives us a value (3 = 0.482. 

7.2 Critical isotherm 

As discussed in the previous section, the critical isotherm at T = T c is the curve marking 
the boundary between single-valued and multi-valued behavior, and correspondingly it has 
a diverging slope for p and s precisely at p = p c . Using the behavior of p on the critical 
isotherm, we can determine the critical exponent 5, defined as 

p-p c ~ \p-p c \ 1/S , for T = T C . (74) 

Now that we have p c in hand, we plot a number of black holes on a log-log plot near the 
critical point with p > p c in figure [8] The points are fit well by a straight line with slope 
giving 5 = 3.03476. We note that the mean field value is 5 = 3. 

7.3 First order axis and susceptibilities 

Susceptibilities in general diverge near a critical point. However, it is not required that all 
possible susceptibilities are divergent. Indeed we find this is the case for a, which is the 
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Figure 8: The rate at which p approaches p c as /i approaches p c on the critical isotherm on 
a log-log plot. The slope gives us a value 5 = 3.03476. 





1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
• 




a = o ; 



-0.003 -0.002 -0.001 0.000 0.001 0.002 0.003 

t 



Figure 9: The specific heat C p near the critical point along a line of constant p, along the 
first-order axis. There is no divergence, giving a = 0. 
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Figure 10: The susceptibility \ 2 near the critical point along a line of constant p. 

power law exponent for the specific heat C p along the axis defined by the first order line: 

C p ~ \T — T c \~ a , along first order axis . (75) 

To avoid the complications of the first-order line itself, we perform this approach from the 
other side, with p, < p c . A calculational advantage is that the line of constant p comes very 
close to the first order axis [51] . so we can simply generate a set of black holes filling out 
that line, and define C p in terms of finite differences of nearest neighbors. 

The result is that C p does not diverge at all along this line, but instead is smooth at 
values near C p ~ 10.5. This corresponds to a vanishing a: 

a = . (76) 



In a sense this is the most robust of all our results, since even a weak divergence looks 



completely different from a lack of divergence; it suggests that the result (76) is exact. 
Moreover, this is again the mean field result; for example in the van der Waals theory of a 
fluid one has C p = 3n/2 with n the number density]^] 

Although we find C p to have no divergence, other quantities do show the expected diver- 
gences. The final thermodynamic exponent, 7, is defined by the approach of \2 along the 
same axis, 

X2 ~ \T — T C \~ J , along first order axis . (77) 



12 Certain systems contain a discontinuity or logarithmic divergence for C p , and are also grouped as a — 0; 
ours is truly smooth, as with the van der Waals field. 
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Figure 11: The baryon susceptibility X2 compared to t = (T — T c )/T c as the critical point is 
approached on a log-log plot. The slope gives us a value 7 = 0.942. 



This exponent is the most difficult to calculate. The quantity involves a derivative in the 
/i-direction, so to calculate the finite difference we must obtain pairs of points with the 
same T to within a very small tolerance AT, which are separated by a larger but still small 
amount A/i; we then need a sequence of pairs of such points moving along the first-order 
axis, a direction unrelated to the derivative. To find black holes near the axis, we again 
imposed constant p. We then looked for pairs of points with Ap < 0.001 and kept those 
with AT/ A// < 0.002. 



The result is shown in figure [T0J The black holes for T — T c > 0.004 form a smooth, 
single-valued curve, but near the critical point numerical errors grow larger and the curve is 
no longer single- valued. Fitting just the single-valued region, we produce the log-log plot in 



figure 11. The resulting exponent is 7 = 0.942, while the mean field value of 7 is ^mf = 1- 



Once again our result is consistent with mean field. 
7.4 Summary and Scaling 

In conclusion, we have measured the critical exponents a, (3, 7 and S and found them to be 
consistent with the mean field values. Since the mean field values are themselves consistent 



with scaling (26), it is clear that our results pass this self-consistency check as well. 



As an idea of the size of the errors in our measurements, we can choose two exponents, 
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use the scaling laws (26) to calculate predictions for the other two, and compare these 
predictions to our actual results. It is natural to use a = as an input, since we obtain it 
as an apparently exact result. Also inputting 5 = 3.035 gives us the results: 



Exponent 


a 


P 


7 


5 


Calculated value 





0.482 


0.942 


3.035 


Scaling prediction from a & S 





0.496 


1.009 


3.035 


% Diff. 




3% 


7% 





Thus our deviations from scaling are in the 3% — 7% range, giving an idea of the size of the 
errors in our method. 

7.5 Other models 

We also carried through an analysis of the model with the same potential ^ and the gauge 
kinetic function 

f(4>) = e-+. (78) 

The phase diagram obtained from this model shares all relevant properties to the one pre- 
sented here: a first-order line ending in a critical point, and critical exponents consistent with 
mean field. We omit the details since they are virtually identical to those just discussed. 
This exponential model has substantially poorer fit to lattice data at fi = 0. 



8 Conclusions 

A number of techniques have been employed to predict the location of the QCD critical 
point. These include lattice calculations that attempt to circumvent the problems of finite 
/i by a number of different means, including taking a Taylor series expansion around fi = 0, 
reweighting the contributions to the path integral, or analytically continuing from imaginary 
chemical potential. There are also calculations in a variety of Nambu-Jona-Lasinio models, 
along with a number of other methods. 



In figure 12, we show the location of a number of different calculations of the location of 
the critical point, along with our result, presented as BH 10 . A key to the various abbreviations 
and references is given in the table; for more information see [281 [32]. The variation in prior 
results is considerable, and our result lies within the parameter space defined by the others. 
Also included in this plot is an estimate for the chemical freeze-out line [5]. 
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Figure 12: Our result for the location of the critical point (BH 10 ) compared to other calcu- 
lations in the literature, along with the chemical freeze-out curve. 



Label 


Method 


Reference 


HB 


Hadronic Bootstrap 


[33] 


LTE 


Lattice Taylor Expansion 


M 


LR 


Lattice Reweighting 


[35], [36] 


RM 


Random Matrix 




NJL 


Nambu-Jona-Lasinio 


[38], [39], [10] 


CJT 


Effective Potential 


m 


LSM 


Linear Sigma Model 


m 


CO 


Composite Operator 





In general, as heavy ion collisions gain center-of-mass energy, the produced medium is charac- 
terized by higher T and lower fi. This leads to several issues with the possibility of exploring 
the region near our prediction, where \x is relatively large. First, assuming the heavy ion 
collisions attain thermodynamic equilibrium, the value of T that can be reached may be 
too small by the time one has reached sufficiently large fx. Moreover, at some large fi the 
center-of-mass energy will become too low to actually thermalize the colliding ions, making 
a thermodynamic interpretation no longer appropriate. Collisions at RHIC have a mini- 
mum energy of 5-7 GeV [5], too high to reach our value of /x c ; LHC is even worse. A more 
promising possibility is the CBM experiment at the future accelerator center FAIR at GSI, 
a fixed-target experiment whose intended region of exploration includes the location of our 
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critical point [BJ. 

That said, as we have emphasized our result should be regarded primarily as a proof 
of principle: it is possible to extract QCD-like phase diagrams from relatively simple holo- 
graphic duals. This result has substantial promise, precisely because finite chemical potential 
calculations are so difficult on the lattice; in gravity duals, finite chemical potential simply 
involves introducing a new field, and possesses no additional qualitative complexity. Of 
course holographic duals introduce other complications — large N and the fact that they are 
not precisely QCD, only QCD-like — that must in turn be dealt with. 

We have discovered a first-order line and a critical endpoint with mean-field critical 
exponents. A number of ways to generalize these results are evident. Most obviously, one 
can study the fluctuations around our classical backgrounds, thereby learning about spectra, 
transport and the true free energy function. An equally obvious, though potentially difficult, 
further step is to add 1/N corrections to the geometries, hopefully moving the critical point 
away from mean field. One can also consider studying a larger theory. In particular, the 
introduction of chiral symmetry in addition to baryon number is straightforward in principle, 
if more intricate in practice. By enlarging the theory one could also hope to study the color 
superconducting phases at large \i. We hope to examine these issues in the future. 
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